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ABSTRACT 


Reliability-based  design  optimization  (RBDO)  seeks  the  best  design  for  a  structural 
system  under  uncertainty.  Typically,  uncertainty  arises  from  random  loads  such  as  wind 
pressure  and  random  material  properties  such  as  yield  stress.  A  reliable  design  must 
account  for  uncertainties  to  ensure  safety. 

Various  methods  have  been  proposed  to  solve  the  nonlinear  optimization  models 
that  RBDO  uses.  However,  these  methods  are  theoretically  and  computationally 
troublesome  as  they  involve  constraints  on  failure  probability,  and  failure  probability  is 
difficult  to  handle  in  optimization  algorithms.  This  thesis  considers  an  alternative 
approach  to  RBDO  that  uses  the  “buffered  failure  probability,”  and  develops  four  new 
solution  algorithms  based  on  sample-average  approximations.  Buffered  failure 
probability  is  more  conservative  than  failure  probability  and  it  is  much  easier  to  handle  in 
optimization  algorithms. 

We  test  the  algorithms  on  six  engineering-design  examples  from  the  literature. 
The  examples  range  from  simple  systems  with  two  design  variables  to  complicated  ones 
with  ten.  Results  show  that  the  new  algorithms  may  reduce  solution  time  by  an  average 
factor  of  560  compared  to  an  existing  algorithm.  Furthermore,  they  can  handle  problem 
instances  with  two  orders  of  magnitude  larger  sample  sizes,  which  may  be  important  for 
reasons  of  accuracy. 
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EXECUTIVE  SUMMARY 


Reliability-based  design  optimization  (RBDO)  deals  with  how  to  best  design  a 
structural  system  under  uncertainty  while  considering  reliability  constraints.  Uncertainty 
present  in  engineering  systems  typically  arises  from  random  loads  such  as  wind  pressure 
and  random  material  properties  such  as  buckling  stress.  A  reliable  design  must  account 
for  these  uncertainties  in  order  to  ensure  safety. 

Various  methods  have  been  proposed  to  solve  the  nonlinear  optimization  models 
that  RBDO  uses.  However,  these  methods  are  theoretically  and  computationally 
troublesome  as  they  involve  constraints  on  failure  probability,  and  failure  probability  is 
difficult  to  handle  in  nonlinear  optimization  algorithms.  The  constraints  based  on  failure 
probability  may  yield  difficult-to-solve  optimization  problems. 

This  thesis  considers  an  alternative  approach  to  RBDO  that  uses  “buffered  failure 
probability,”  and  considers  five  solution  algorithms  based  on  sample-average 
approximations.  The  buffered  failure  probability  approach  is  more  conservative  than  the 
traditional  approach  based  on  the  failure  probability,  meaning  that  a  design  that  satisfies  a 
reliability  constraint  based  on  buffered  failure  probability  is  guaranteed  to  satisfy  one 
based  on  failure  probability.  The  buffered  failure  probability  is  also  much  easier  to  handle 
in  nonlinear  optimization  algorithms. 

The  first  three  algorithms  each  solve  models  that  incorporate  a  single  sample- 
average  approximation  of  the  buffered  failure  probability  constraint.  Algorithm  1  is  a 
well-known  method  based  on  a  model  refonnulation  and  solves  a  single,  large  nonlinear 
program.  Algorithms  2-5  are  new  algorithms  developed  in  this  thesis.  Algorithm  2  is  an 
active-set  implementation  of  Algorithm  1 .  Algorithm  3  uses  exponential  smoothing  of  a 
max-function  to  avoid  the  large-scale  reformulation  of  Algorithm  1.  Algorithm  4 
approximately  solves  a  sequence  of  sample-average  approximations  within  an  adaptive 
sample-adjustment  scheme  that  ensures  the  sample  size  is  gradually  increased  to  infinity. 
Algorithm  5  is  similar  to  Algorithm  4,  but  includes  an  active-set  strategy. 


xv 


We  test  the  algorithms  on  six  engineering-design  examples  from  the  literature. 
The  examples  range  from  simple  systems  with  two  design  variables  to  complicated  ones 
with  ten  design  variables.  Results  from  Algorithm  2  show  an  average  speed-up  in 
solution  time  by  a  factor  of  560  over  Algorithm  1.  Algorithm  3  exhibits  a  speed-up  by  a 
factor  of  3 1  over  Algorithm  1 .  Algorithms  2  and  3  can  handle  problem  instances  with 
large  sample  sizes,  in  fact,  one  and  two  orders  of  magnitude  larger  than  that  of  Algorithm 
1,  respectively.  The  ability  to  handle  large  sample  sizes  is  important,  for  reasons  of 
accuracy.  Algorithms  4  and  5  obtain  high-quality  solutions  in  minutes  without  the  need 
for  a  user  to  specify  a  sample  size,  which  may  be  difficult  in  practice. 

The  results  in  this  thesis  show  that  Algorithms  2-5  have  significantly  improved 
engineers’  ability  to  relatively  quickly  generate  cost  efficient  designs  that  satisfy  a  failure 
probability  constraint. 


xvi 


ACKNOWLEDGMENTS 


I  would  like  to  extend  my  personal  thanks  to  Professor  Johannes  O.  Royset,  for 
his  clear  guidance  and  constant  availability.  With  his  unwavering  patience  this  difficult 
process  has  become  an  exceptional  learning  opportunity.  I  could  not  have  completed  this 
research  without  his  tireless  efforts. 

I  also  wish  to  thank  Professor  R.  Kevin  Wood  for  his  dedicated  efforts  to  ensure 
this  study  to  be  genuine  research  writing.  His  editing  contribution  and  professional 
insight  ensured  the  highest  standard  in  this  thesis. 

I  want  to  thank  my  wife,  Medalet,  for  her  unconditional  love,  precious  support 
and  endless  patience  in  the  midst  of  this  challenge  and  my  daughter,  Elif,  for  enduring 
days  without  dad  so  I  could  finish  this  work. 

Finally,  I  would  like  to  offer  my  special  thanks  to  my  country,  Turkey,  for 
providing  me  the  unique  opportunity  to  attend  the  Operations  Research  program  at  NPS. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

This  thesis  considers  the  reliability-based  design-optimization  problem  (RBDO) 
of  minimizing  the  design  cost  of  a  structural  system  subject  to  a  constraint  on  its 
reliability  and  other  quantities.  We  characterize  a  system’s  reliability  in  terms  of  the 
probability  of  failure  with  respect  to  one  or  more  performance  requirements. 

We  approximate  this  problem  by  replacing  the  failure  probability  constraint  with 
a  buffered  failure  probability  constraint  and  develop  four  algorithms  for  solving  the 
resulting  nonlinear  program. 

B.  MOTIVATION 

Engineers  aim  to  minimize  design  costs  of  structural  systems  such  as  aircraft 
wings,  vehicle  frames,  ship  hulls,  bridges,  and  buildings.  The  minimization  typically  is 
subject  to  one  or  more  constraints  on  system  reliability.  The  parameters  that  describe  the 
shape  and  characteristics  of  the  system  are  referred  to  as  design  variables.  These  variables 
are  controlled  by  the  engineer  and  are  manipulated  so  that  the  cost  of  the  system  is 
minimized  and  the  reliability  is  sufficiently  large.  Therefore,  the  challenge  for  an 
engineer  is  to  reduce  the  design  cost  while  satisfying  requirements  for  system 
performance  and  reliability. 

A  system’s  reliability  accounts  for  uncertain  loads  that  act  on  the  system,  and  the 
uncertain  capacity  of  the  system  to  withstand  these  loads.  There  are  a  variety  of 
uncertainties  that  engineers  need  to  consider.  A  system’s  performance  depends  on  the 
type  and  magnitude  of  loads  and  the  strength  of  the  system,  which  is  related  to  properties 
of  materials  used  in  the  design.  The  loads  and  system  strength  can  be  described  by 
random  (uncontrollable)  variables.  As  an  example,  for  a  building,  wind  pressure  is  an 
uncontrollable  load  that  can  be  modeled  using  random  variables.  The  loads  may  lead  the 
system  to  not  being  able  to  meet  functional  and  safety  requirements.  The  result  can  be 
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loss  of  serviceability,  or  a  complete  destruction  of  the  system.  The  functional 
requirements  are  called  limit  states  of  the  system  and,  when  they  are  exceeded,  we  say 
that  the  system  is  failed. 

Various  methods  have  been  proposed  for  solving  the  nonlinear  programs  arising 
in  RBDO.  However,  these  methods  are  theoretically  and  computationally  troublesome  as 
they  involve  the  failure  probability,  which  may  be  nonsmooth  and  nonconvex  and 
therefore  difficult  to  handle  by  standard  nonlinear  programming  solvers  such  as  SNOPT 
(Gill,  Murray,  &  Saunders,  1998),  LANCELOT  (Conn,  Gould,  &  Toint,  1992),  and 
NLPQL  (Schittkowski,  1985).  This  thesis  explores  another  approach  to  RBDO  based  on 
buffered  failure  probability.  The  buffered  failure  probability  approach  is  more 
conservative  than  the  traditional  approach  based  on  the  failure  probability  (Rockafellar  & 
Royset,  2010).  Thus,  a  design  that  satisfies  the  reliability  constraint  based  on  the  buffered 
failure  probability  also  satisfies  one  based  on  the  failure  probability.  The  buffered  failure 
probability  is  computationally  easier  to  handle  due  to  the  algorithmic  advances  of  this 
thesis.  Rockafellar  and  Royset  (2010)  also  discuss  other  potential  advantages  that 
buffered  failure  probability  has  over  failure  probability  in  this  context.  However,  this 
thesis  shows  that  the  computational  advantages  alone  are  sufficient  to  warrant  a 
preference  for  the  buffered  failure  probability  approach. 

C.  SCOPE  AND  LIMITATIONS 

We  represent  uncertainties  as  random  variables.  In  principle,  the  random  variables 
can  be  either  discrete  or  continuous.  However,  we  only  consider  continuous  random 
variables  in  this  thesis. 

We  assume  that  we  know  the  joint  distribution  of  the  random  variables  and  that 
the  corresponding  cumulative  distribution  function  is  strictly  increasing. 

A  real-life  structural  system  generally  is  composed  of  several  components.  If 
failure  of  any  one  of  these  components  constitutes  failure  for  the  entire  system,  then  we 
call  the  structure  a  series  system.  In  contrast,  parallel  systems  are  those  that  need  the 
failure  of  all  the  components  for  a  system  failure.  In  this  thesis,  we  focus  on  series 
systems  with,  consequently,  one  failure  probability  constraint. 
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D. 


LITERATURE  REVIEW 


The  classical  methods  for  solving  the  RBDO  perform  a  double-loop  approach:  an 
outer  optimization  loop  and  an  inner  reliability  assessment  loop;  see  Du  and  Chen  (2004). 
The  reliability  assessment  is  performed  by  two  different  methods:  Reliability  Index 
Approach  (RIA)  (Lee,  Yang,  &  Ruy,  2002)  or  Performance  Measure  Approach  (PMA) 
(Tu,  Choi,  &  Park,  1999).  Both  of  these  methods  approximate  the  failure  probability  by 
surrogates  of  unknown  accuracy  and  hence  cannot  guarantee  the  convergence  to  globally, 
locally,  or  stationary  solutions  of  RBDO  problems.  A  single-loop  RBDO  approach  is 
proposed  by  Liang,  Mourelatos,  and  Tu  (2008).  This  method  utilizes  the  fact  that  a 
surrogate  for  the  failure  probability  can  be  evaluated  by  solving  a  nonlinear  program. 
Hence,  the  RBDO  problem  with  a  failure  probability  constraint  can  be  approximated  by 
an  optimization  model  with  equilibrium  constraint.  However,  the  resulting  model  may  be 
difficult  to  solve  due  to  nonconvexity,  nonsmoothness,  and/or  lack  of  a  constraint 
qualification.  Moreover,  as  the  accuracy  of  the  surrogate  for  the  failure  probability  is 
unknown,  the  computed  design  may  be  nonoptimal  and  violate  failure  probability 
constraints.  We  refer  to  Rockafellar  and  Royset  (2010)  and  Royset,  Der  Kiureghian,  and 
Polak  (2006)  for  a  more  comprehensive  literature  review  and  difficulties  associated  with 
current  approaches. 
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II.  BACKGROUND  AND  PROBLEM  DEFINITION 


A.  LIMIT-STATE  FUNCTIONS  AND  FAILURE  PROBABILITY 

A  system  response  relative  to  a  functional  requirement,  viewed  as  a  function  of 
the  design  and  random  variables,  is  referred  to  as  a  limit-state  function  and  denoted  by 
g(x,V).  Here  x  =  (xl,x2,...,xn)'  is  a  vector  of  design  variables  (where  prime  '  denotes 

the  transpose  of  a  vector)  and  V  =  (VvV2,...,Vm)'  is  a  vector  of  random  variables.  The 
joint  probability  distribution  of  the  random  variables  is  known.  We  denote  realizations  of 
V  by  v.  For  a  given  x,  g(x,v)  represents  the  perfonnance  of  the  system  under 
realization  v  of  V.  An  unsatisfactory  performance,  i.e.,  “failure,”  occurs  when 
g(x,v)  >  0.  If  g(x,v)  <  0  then  the  system  perfonnance  is  acceptable. 

In  complex  systems,  there  may  be  multiple  limit-state  functions;  for  example,  see 
Example  5  in  the  Appendix,  taken  from  Rao  (2009,  pp.  472-473).  We  denote  these 
functions  as  gk(x,v),  k  ef  ,  where  K  =  {1,2 ,...,/«}.  We  also  define 

g(x,v)  =  max{gi(x,v)}.  (1) 

keK 

We  can  characterize  the  reliability  of  the  system  by  its  failure  probability  defined 
by 

p(x)  =  Prob[g(x,V)>  0].  (2) 

The  following  simplified  example  illustrates  a  design  process  based  on  the  failure 
probability.  (Note:  For  clarity,  this  example  uses  discrete  x,  but  we  note  that  this  thesis 
develops  solution  methods  only  for  continuous  x. ) 

Example  1.  A  TOW  missile  is  a  heavy  anti-tank  missile.  One  component  of  the 
missile’s  launcher  is  an  optical  system.  Suppose  that  two  different  optical  systems,  1  and 
2,  are  available  for  incorporation  in  a  final  design:  setting  xt  =  1  means  system  i  is 

selected  and  x;  =0,  otherwise.  The  decision  maker  plans  to  procure  and  incorporate  the 

system  with  lower  failure  probability.  Systems  1  and  2  have  an  operational  capability 
down  to  -28  °F  and  -30  °F,  respectively.  The  procurement  cost  of  system  2  is  higher 
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than  system  1 .  The  missile  will  be  operated  under  severe  conditions  in  a  certain  mission 
area.  Let  V  represent  the  random  temperature  in  Fahrenheit  degrees,  which  is  known  to 
be  normally  distributed  with  mean  -20  °F  and  standard  deviation  3  °F  in  that  area. 
Engineers  have  characterized  the  reliability  of  the  system  by  the  following  limit-state 
function 

g(x,v)  = -0.9v-28xj -30x2,  (3) 

where  xl  +  x2  =  1  and  xnx2  e  {0,1}.  The  coefficient  0.9  is  a  scaling  factor  and  represents 
the  impact  of  the  current  temperature  on  the  system.  If  g(x,v )  is  positive  for  a  given 
temperature  v,  the  system  fails.  We  compute  that  designs  1  and  2  have  failure 
probabilities  of  1.06x10  4  and  4.406xl0~6,  respectively.  Thus,  system  2  will  be 
procured. 

For  a  given  design  x ,  the  computation  of  failure  probability  requires  the 
evaluation  of  a  high-dimensional  integral.  That  is, 

p(x)  =  J...J  I  (g(x,v))fv(V)dvv..dvm,  (4) 

where  fv(V)  is  the  joint  probability  density  function  for  the  random  vector  V,  and 
I(z)  = 1  if  z  >  0,  and  /(z)  =  0  otherwise. 

B.  DESIGN  OPTIMIZATION  OF  A  SYSTEM  SUBJECT  TO  A  FAILURE 

PROBABILITY  CONSTRAINT 

Engineers  often  seek  to  solve  the  design-optimization  problem  (Rockafellar  & 
Roy  set,  2010) 

P :  min  fix) 

A 

s.t.  p(x)<\-a  (5) 

X  G  X, 

where  f(x)  is  a  detenninistic  and  continuously  differentiable  cost  function  for  the 
system,  X  is  a  continuous  region  of  allowable  designs,  and  a  is  a  desired  reliability 
level  in  (0,1].  We  also  assume  that  X  is  convex. 
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Problem  P  represents  the  current  approach  to  RBDO.  Typically,  P  is  difficult  to 
solve  even  approximately  because  the  computation  of  the  failure  probability  in  (4)  is 
computationally  challenging.  Furthermore,  the  integrand  in  (4)  is  non-smooth,  i.e.,  is  not 
differentiable  at  every  point,  and  this  may  cause  difficulties  during  optimization. 
Moreover,  p(x)  is  generally  nonconvex.  With  lack  of  convexity,  a  solution  algorithm 
may  yield  only  a  locally  optimal  solution. 

Solving  P  is  also  difficult  because  computation  of  the  failure  probability  requires 
the  evaluation  of  the  high-dimensional  integral  in  (4).  For  a  real-world  system,  the  limit- 
state  functions  are  highly  complex  and  involve  many  different  random  variables.  (See 
Example  6  in  the  Appendix,  taken  from  Samson,  Thoomu,  Fadel,  &  Reneke  2009). 
Consequently,  failure  reliability  can  only  be  estimated:  the  standard  approach  uses  Monte 
Carlo  simulation  (MCS)  (Melchers,  1999;  Choi,  Grandhi,  &  Canfield,  2007). 

MCS  is  a  sampling  method  that  estimates  expectation  and  probability.  When 
estimating  the  probability  of  failure,  MCS  computes  the  following  estimate  of  the 
probability  of  failure  in  (4): 

1  N 

P(x)  =  —YiI(g(x,vJ)),  (6) 

A  j=i 

where  N  is  the  number  of  random,  sampled,  vector  realizations  vj .  The  function  p(x)  is 
an  unbiased  estimator  of  p(x)  (Rubinstein  &  Kroese,  2008).  The  standard  deviation  of 
p(x)  is  inversely  proportional  to  the  square  root  of  the  sample  size  N.  Hence,  the 
accuracy  of  failure  probability  estimates  is  poor  for  small  samples,  especially  for  small 
failure  probabilities.  Therefore,  replacing  the  failure  probability  by  its  estimator  in  (6) 
computed  with  a  large  N  leads  to  long  solution  times  for  sampled  approximations  to  P. 
Furthermore,  (6)  is  not  differentiable  because  of  the  indicator  function,  and  this  may 
cause  convergence  difficulties  for  standard  nonlinear  optimization  algorithms.  We  refer 
to  Rockafellar  and  Royset  (2010)  for  more  detailed  explanations  of  difficulties  associated 
with  p(x). 
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c. 


BUFFERED  FAILURE  PROBABILITY 


As  discussed  above,  the  failure  probability  has  several  undesirable  properties  that 
may  cause  difficulties  for  a  standard  nonlinear  optimization  algorithm.  We  now  introduce 
an  alternative  approach  to  the  design-optimization  problem  using  another  measure  of 
failure,  the  buffered  failure  probability  (Rockafellar  &  Royset,  2010).  The  buffered 
failure  probability  approach  is  based  on  the  conditional  value-at-risk  (Rockafellar  & 
Uryasev,  2000;  Rockafellar  &  Uryasev,  2002),  which  is  used  in  financial  engineering  to 
compute  optimal  investment  portfolios.  We  next  define  the  buffered  failure  probability. 

Suppose  we  wish  to  solve  P  with  a  failure  probability  level  l- a.  Instead  of 
imposing  the  constraint  p(x )  <  1  -  a,  we  introduce  p(x)  <  1  -  a  as  an  alternative, 
where  p(x )  denotes  the  buffered  failure  probability  (Rockafellar  &  Royset,  2010).  A 
design  that  satisfies  the  buffered  failure  probability  constraint  also  satisfies  the  failure 
probability  requirements,  as  we  describe  later,  and  the  constraint  p(x)  <  1  -  a  is  easier 
to  handle  computationally.  The  following  description  follows  Rockafellar  and  Royset 
(2010)  closely. 


Let  the  cumulative  distribution  function  of  g(x,  V)  be  denoted  by 

Ffy)  =  Prob(g(x,V)<y),  (7) 

which  we  assume  to  be  continuous  and  strictly  increasing.  We  next  define  the  a-  quantile 
function,  or  simply  a-  quantile,  which  is  denoted  as  qa(x).  For  a  probability  level  a, 


In  addition,  we  define  the  a- super  quantile  function,  or  simply  a-  superquantile,  by 

fa  O)  =  E  [six,  V)  |  g(x,  V )  >  qa  (x)] . 

Integration  is  also  useful  in  finding  the  superquantile.  That  is, 

<?«(*)  = 


1  -a 


\qp(x)dp. 


(8) 

(9) 


(10) 


The  superquantile  can  be  interpreted  as  the  value  of  g(x,  V )  that  splits  the  area  under  the 

probability  density  function  in  the  interval  [ga(x),co)  into  two  balancing  parts.  That  is,  it 

represents  the  average  value  of  the  limit-state  function,  conditioned  on  g(x,  V )  being  no 
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less  than  the  a-  quantile.  The  superquantile  is  identical  to  the  conditional  value-at-risk, 
but  we  follow  Rockafellar  and  Royset  (2010)  here  and  use  the  term  superquantile.  When 
the  superquantile  equals  zero  at  a  probability  level  a ,  then  the  buffered  failure 
probability  p(x )  is  equal  to  l  -  a.  Hence,  the  buffered  failure  probability  may  be 
defined  as 


p(x)  =  Prob  [g(x,  V )  >  qaQ  (x)] , 


(11) 


where  a0  is  such  as  <7  (x)  =  0  (Rockafellar  &  Royset,  2010). 


Figure  1  illustrates  the  probability  density  function  for  an  example  limit-state 
function  at  a  fixed  x  =  x,  where  a0  has  been  identified  such  that  q  (x)  =  0.  qaQ  (x)  splits 

the  area  under  the  probability  density  function  to  the  right  of  -2  into  two  balancing  parts. 
Hence,  the  quantile  qao  (x)  =  -2 .  We  also  see  in  Figure  2,  which  represents  the 

cumulative  distribution  function  for  the  same  limit-state  function,  that  a0  therefore  must 
equal  0.85.  Thus,  we  calculate  the  buffered  failure  probability  for  this  example  as 
1-0. 85, which  is  0.15. 


g(x,v) 

Figure  1.  Example  Probability  Density  Function  (pdf)  for  g(x,V )  when  That 

Function  is  Normally  Distributed 
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F(g(Xv)) 


Figure  2.  Example  Cumulative  Distribution  Function  (cdf)  for  g(x,  V)  when  That 

Function  is  Normally  Distributed 

In  general,  for  any  a  value  and  any  xeX,  qa(x)<qa(x)  and  thus 
p(x)  <  p(x).  Flence,  the  buffered  failure  probability  is  a  conservative  estimate  of  the 
failure  probability.  The  relationship  p(x)  <  p(x)  is  easily  obtained  from  the  following 
argument.  Suppose  x  =  x  is  fixed,  and  a0  is  chosen  so  that  qaQ  (x)  =  0.  Then,  by  the 
definition  of  the  superquantile,  it  follows  that  qBo  (x)  <  0.  Since  q  (x)  =  0,  the  buffered 
failure  probability  p(x)  equals  to  1  -a0.  It  follows  from  the  definition  of  the  failure 
probability  that  if  qa<(x)  =  0  for  some  av  then  p(x)  =  \-av  Since  qa()  (x)  <  qfA  (x),  it 
must  be  true  that  a0  <ay  Thus,  1  -a0  is  no  smaller  than  1  -ar  Since  the  arguments 
above  hold  for  any  x  eX  we  have  the  result  that  p(x)  <  p(x)  for  all  xel.  Figure  2 
illustrates  this  situation  with  ax  =  0.96,  a0  =  0.85,  q  (x)  =  0,  and  qao  (x)  =  -2. 

The  computation  of  the  buffered  failure  probability  appears  cumbersome.  As  we 
see  below,  however,  we  never  directly  use  this  definition  in  computations  and  instead  use 
alternative  expressions  easily  incorporated  in  optimization  models. 
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The  buffered  failure  probability  constraint  p(x)  <  1  -  a  is  satisfied  if 


qa(x)<0. 

Rockafellar  and  Uryasev  (2002)  show  that  for  any  rel 

qa(x)  =  m\n<f>a(z(j,x), 

zo 

where  z0  is  an  auxiliary  design  variable  and 

4  Oo  ,x)  =  z0+  — E  [max  { 0,  g(x,V)-z0}]. 

1  -a  L  J 

Therefore,  a  design  rel  and  a  value  for  z0  that  satisfy 

<f>a(z0,x)<0 

also  satisfy  p(x)  <  1  -  a. 


(12) 

(13) 


(14) 


(15) 


Although  p(x)  generally  overestimates  p{x),  the  numerical  examples  that  we 
discuss  in  Chapter  IV  show  that  the  difference  between  the  two  probabilities  need  not  be 
great. 

We  can  formulate  a  restriction  of  P,  by  replacing  the  reliability  constraint  in  P 
with  (15).  Thus,  a  restriction  to  P  using  buffered  failure  probability  takes  the  form 
BP  :  min  / (x) 

W  z0 

s.t.  Z0+— !—  £fmax{0,g(x,F)-z0}]<0  (16) 

1  —  a 

x  €  X. 

With  BP  being  a  restriction  of  P,  the  solution  of  BP  is  also  feasible  in  P,  but  may  not  be 
optimal  in  P.  However,  the  constraint  (16)  is  easier  to  handle  than  (5)  as  it  is  convex 
when  each  gk{x,v),  keK  is  convex  in  x  for  all  v  (Rockafellar  &  Royset,  2010). 

Moreover,  in  contrast  to  the  nonsmoothness  of  the  indicator  function  in  (4),  the 
nonsmoothnesss  in  (16)  is  easily  handled  (as  described  in  Chapter  III).  Thus,  computation 
of  optimal  values  is  easier  in  BP  than  in  P.  We  also  note  that  the  buffered  failure 
probability  accounts  for  the  tail  of  the  distribution  of  g(x,  V)  in  a  different  way  than  does 
failure  probability.  We  refer  to  Rockafellar  and  Royset  (2010)  for  a  discussion  of  why 
this  might  be  an  advantage  for  RBDO. 
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We  next  illustrate  the  failure  probability  and  buffered  failure  probability 
approaches  on  a  stochastic,  continuous  knapsack  problem  with  random  knapsack 
capacity.  Let  c  and  w  be  deterministic  vectors  of  per-unit  item  values  and  item  weights, 
respectively,  let  or  be  a  user-defined  reliability  level,  let  f  be  a  random  variable 
describing  the  total  capacity  of  the  knapsack  in  terms  of  weight,  and  let  x  be  a  vector  of 
decision  variables  that  chooses  the  amount  of  each  item  to  placed  in  the  knapsack.  We 
assume  that  V  is  normally  distributed  with  mean  //  and  standard  deviation  a .  The 
knapsack  problem  then  takes  the  fonn: 

KP :  max  c'x 

x 

s.t.  Prob(w'x  <V)>  a  (17) 


x  >  0. 


In  our  notation,  this  problem  instance  has  only  one  limit-state  function: 

g(x,V)  =  w'x-V.  (18) 

Since  limit-state  function  values  greater  than  zero  represent  failure,  the  following 
equation  becomes  the  reliability  constraint: 

Prob(w'x-V  >  0)  < \-a,  (19) 

where  w'x-V  is  normally  distributed  with  mean  w'x  -  //  and  standard  deviation  a . 
The  reliability  constraint  in  ( 1 9)  takes  the  following  form: 

(0-(wx-V)) 


1-0 


<  1  -a, 


(20) 


J 


v  cr 

where  O(-)  represents  the  cumulative  distribution  function  of  the  standard  normal 
distribution.  Thus,  the  problem  KP  is  equivalent  to 


KP  : 


mm  —  c  x 
x 


s.t.  w'x-ju  +  cj  O  l(a)  <  0 


(21) 


x  >  0. 

Next,  we  consider  the  buffered  failure  probability  approach.  In  this  case  BP  takes  the 
form 
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KBP: 


min  -cx 
x 


s.t.  qa(x)<  0 


(22) 


x  >  0. 

Since  g(x,V )  is  nonnally  distributed  with  mean  w'x-  ju  and  standard  deviation  cr,  the 
superquantile  qa{x)  is  easily  computed  (Truncated  Normal  Distribution,  Wikipedia,  n.d.) 
and  we  obtain 


qa(x)  =  n  + 


(23) 


where  (/>(■)  is  the  standard  normal  probability  density  function.  Therefore,  KBP  takes  the 
following  equivalent  fonn: 

KBP':  min  -cx 

x 


s.t. 


w'x  -ju  + 


I  -  a 


(24) 


x  >  0. 

Thus,  for  this  knapsack  problem,  the  difference  between  the  reliability  constraints  of  the 
failure  and  buffered  failure  probabilities  corresponds  to  the  difference  between  (21)  and 
(24).  Table  1  gives  the  near-optimal  solutions  of  KP'  and  KBP'  when  V  is  normally 
distributed  with  mean  3.5  and  standard  deviation  0.1,  c'  =  (2,l),  w'  =  (1.1, 2.1) ,  and 
a  =  0.99  .  From  Table  1  we  see  that  the  buffered  failure  probability  approach  results  in  a 
conservative  design,  with  a  failure  probability  that  is  about  40%  lower  than  required. 


Problem 

Xi 

x2 

f 

-cx 

p(x) 

KP' 

1.06124 

1.00000 

-3.12248 

0.0100 

KBP' 

1.03043 

1.00000 

-3.06087 

0.0039 

Table  1  Comparison  of  Solutions  Using  Failure  and  Buffered  Failure  Probabilities 
for  a  Knapsack  Problem  with  V  ~  N(3. 5,0. 12),  c'  =  (2,l),  w'  =  (1.1, 2.1),  and 

a  =  0.99 
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III.  ALGORITHMS  FOR  RBDO  BASED  ON  BUFFERED  FAILURE 

PROBABILITY 


This  chapter  examines  five  algorithms  that  exploit  the  properties  of  buffered 
failure  probability  to  compute  an  approximated  solution  of  BP  in  a  relatively  short  time 
and  with  high  accuracy.  Since  BP  is  a  restriction  of  P,  each  solution  also  represents  an 
approximate  solution  of  P.  Algorithm  1  is  a  well-known  method  based  on  the  solution  of 
a  reformulation  of  BP.  Algorithms  2-5  are  new  algorithms  developed  in  this  thesis. 
Algorithm  2  is  an  active-set  strategy  implementation  of  Algorithm  1.  Algorithm  3 
implements  exponential  smoothing  of  the  max-function  in  (16)  and  proceeds  to  solve  a 
smoothed  problem.  Algorithms  1-3  approximate  the  expectation  A[max  JO,  g(x,  v)  -  z0 }  ] 
in  (16)  by  its  sample  average  for  a  fixed  sample  size  N ,  i.e., 

1  N 

—  ^max{°,g(x,vy)-Zo}  (25) 

A  j= i 

with  random  vector  realizations  vj ,  j  =  1,2,...,  A  .  Thus,  the  reliability  constraint  in  (16) 
takes  the  following  form: 

+  N^_  —  X  max  {°>  S(x,  vj)  -  z0 }  <  0.  (26) 

Algorithm  4  implements  an  adaptive  sample-adjustment  scheme  that  ensures  the 
sample  size  is  gradually  increased  to  infinity.  Iterates  generated  by  Algorithm  4  are 
guaranteed  to  converge  to  a  stationary  point  of  BP  (Royset,  2010b).  Algorithm  5  is 
similar  to  Algorithm  4  but  includes  an  active-set  strategy.  We  next  discuss  each  algorithm 
in  turn. 

A.  ALGORITHMS  WITH  FIXED  SAMPLE  SIZE 

1.  Algorithm  1:  Reformulation  of  BP  Using  Sampling  and  Auxiliary 
Variables 

The  constraint  (26)  is  non-smooth  and  is  therefore  not  tractable  by  standard 
nonlinear  optimization  algorithms.  Algorithm  1,  an  existing  algorithm,  uses  a 
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reformulation  to  overcome  this  difficulty.  This  refonnulation  was  first  proposed  in 
Rockafellar  and  Royset  (2010). 

We  introduce  auxiliary  variables  z,,  j  =  ,  and  find  that  for  any  j, 

max|0,g(x,v7)-z0}  <  z.  (27) 

is  equivalent  to 

0  <  z .  (28) 

and 

g(x,vJ)-z0<Zj.  (29) 

Hence,  we  can  use  (28)  and  (29)  to  reformulate  the  “max”  part  of  the  constraint  in  (26). 
This  results  in  the  following  transcription: 


BP N  :  min  f(x) 

x,z 

1  N 

s.t.  z  + - y  z  <  o 

N(l-a)jl  J 

(30) 

gk(x,VJ)-Z0<Zj, 

V /  eJ  \/k  ef 

(31) 

JNI 

IV 

o 

VyeJ 

(32) 

xeX, 

where  z  =  (z0,zl,...,zN)'  and  J  =  (1,2, ...,7V} .  For  large  N,  BPjv  is  approximately 
equivalent  to  BP  (Rockafellar  &  Royset,  2010).  However,  BPy  necessitates  introducing  a 
new  auxiliary  design  variable,  i.e.,  z; ,  for  every  realization  of  the  random  vector:  This 

typically  results  in  large-scale  problems.  Given  BPat,  our  first  algorithm  takes  the 
following  simple  form: 

Algorithm  1:  Apply  a  standard  nonlinear  solver  to  BP^  . 

We  use  the  SNOPT  (Gill  et  al.  1998)  solver  here  and  below  in  the  following  algorithms 
in  computational  experiments. 
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2.  Algorithm  2:  External  Active-Set  Strategy 

BP#  becomes  a  large-scale  problem  as  N  increases.  Algorithm  2  aims  to 
overcome  this  difficulty  by  using  an  active-set  strategy  proposed  in  Chung,  Polak,  and 
Sastry  (2010). 

In  this  approach,  we  do  not  let  the  standard  nonlinear  solver  consider  constraints 
and  variables  assumed  to  be  “unimportant”  at  an  optimal  solution.  The  following  is  an 
adaptation  of  the  algorithm  proposed  in  Chung  et  al.  (2010)  to  BP^.  Let  y  =  (x.z)  and 
8  >  0 .  Then,  we  let 

ffl i(y)  =  gk(x,vi)-za-zJ ,  V jeJ  \/k  e  K  (33) 


Z(y)=  max  a>i(y) 

j<=J,keK 


(34) 


Z(y)+  ~  max{0,  z(y)} 


(35) 


we{y)  =  (O', k)  e  J x K  |  o)l(y)  >  %(y)+  - ,  (36) 

where  w£(y )  represents  “active”  constraints  at  y  .  In  each  iteration  of  this  algorithm  we 

limit  the  number  of  iterations  of  the  solver.  Let  denote  the  solution  from  n 

iterations  of  the  solver,  starting  from  y ,  applied  to  the  following  problem: 

BP N(W):  min  f{x) 

X,Z 

1  n  (37) 

s.t.  z„  H - >  z  <  0 

J 

gk(x,vJ)-z0<Zj,  (j,k)  e  W  (38) 

z7  >  0,  V/  (39) 

xel. 

We  note  that  for  j  eJ  such  that  ( j,k )  g  W  for  any  k  e  K  ,  z.  in  BP N(W)  can  be  set  to 
zero. 

Algorithm  2: 

Data:  yQ  initial  guess  for  variable  values,  s>  0 ,  o  >  0 ,  and  n  >  0  an  integer. 
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Step  0:  Set  i  —  0  and  initial  working  set  W.  =  we(y.) . 

Step  1:  Compute  yM  =  A^(y). 

Step  2:  Compute  %(yM) . 

if  z(yi+i) - o  and  \z(yM)-z(yi)\^v 

STOP, 

else  compute  ws(yi+]), 

WM=W,v{wt(yM)}, 

Replace  i  by  i  + 1  and  go  to  Step  1 . 

In  Algorithm  2,  we  hope  to  reduce  the  solution  time  by  considering  only  the 
active  constraints  in  the  optimization. 


3.  Algorithm  3:  Exponential  Smoothing  of  Max-Function 


We  next  introduce  exponential  smoothing  for  BPW.  Let  L=  {0,l,2,...,m}  and 


g0(x,v)  =  0  for  all  x  and  v.  Then,  the  following  equations  represent  the  exponential 
smoothing  of  the  max-function  in  (26)  (Polak,  Womersley,  &  Yin,  2008): 

nJk(x,z0)  =  gk(x,vJ)-z0,  Vy  e  J  VkeL  (40) 

CJ  (x,  z0)  =  max  ijl  (x,  z0),  VyeJ  (41) 

k<=L 

and 


f 


$p(z  zo)  =  CJ  (x,z0)  +  -log 


m  p 

z« 

k= 0 


riJk(x,z())-Cj  (x,z«) 


V/'eJ 


(42) 


V  J 

where  p  >  0  is  a  smoothing  parameter.  Thus,  replacing  the  reliability  constraint  in  (26) 
by  (42)  we  obtain  the  following  approximation  to  BPa?: 
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(43) 


BIV  min  f{x) 

l  N  ■ 

z0  + - V  SJn(x,zQ)<  0 

0  N(l-ct)P  p  ° 

x  e  X. 

If  all  the  limit-state  functions  are  continuously  differentiable,  which  we  assume  they  are 
in  this  thesis,  then  SJp(x,z0),  j  ei,  are  also  continuously  differentiable. 

We  improve  the  quality  of  the  smoothing  approximation  as  we  increase  the  value 
of  p  .  That  is,  the  error  in  the  constraint  due  to  smoothing  vanishes,  as  p  — >  oo  (Polak, 
Womersley,  &  Yin,  2008).  Our  third  algorithm  then  takes  the  following  form: 

Algorithm  3:  Apply  a  standard  nonlinear  solver  to  BP^p. 


B.  ALGORITHMS  WITH  INCREASING  SAMPLE  SIZE 


We  next  introduce  another  approach,  which  adaptively  increases  the  sample  size 
during  the  optimization  process.  This  approach  intends  to  avoid  using  unnecessarily  large 
sample  sizes.  Suppose  that 

X  =  {x\fJ(x)<0,  j  =  1,2,..., <7},  (44) 

where  f1  (x)  represents  nonlinear  and  non-negativity  constraints  and 


1  N  ■ 

$Np(X,  Z0)  =  Z0  +  NQ_a)  Z  SP  Z«) 

represents  the  reliability  constraint  defined  in  (43).  Moreover,  let 

^(x,z0)=max{^(x,z0),  max  fj(x)  J, 
Ynp  (A  z0  )+  =  max  { 0,  y/Np  (x,  z0 )} , 
and 

0Np(x,z  0)  =  min{ma  x{-i//Np(x,z0)++Vf(x)'h, 


max[^  (Y  z0 )  -  Ynp  (*>  Zo)+  +  V^Np  (x,  zJh  +  (x,  z0 )'£ 


(45) 

(46) 

(47) 


max 

j= K-,q 


{/'  W  -  V„  (*,  z, ).  +  Vp  (x)’K)  ] }+ i  ||*f  +  Jp  } 


(48) 
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We  note  that  the  calculation  of  0N ■  (x,z0)  requires  the  solution  of  a  convex  quadratic 

program  and  can  therefore  be  carried  out  in  finite  time  (Royset,  2010a).  We  use  these 
tenns  in  the  definition  of  Algorithms  4  and  5  below. 

1.  Algorithm  4:  Adaptive  Sample-Size  Algorithm 

Algorithm  4  is  similar  to  one  described  in  Royset  (2010a),  but  includes  the  use  of 
smoothing;  see  Royset  (2010b).  Let  A"Np(x,z0)  denote  the  solution  from  n  iterations, 

starting  from  x  and  z0 ,  of  a  standard  nonlinear  solver  applied  to  BP  V/). 

Algorithm  4: 

Data:  /V0  initial  sample  size,  p0>  0  initial  smoothing  parameter,  s  e  (0.1]  sample-size 
increase  factor,  xQ  initial  design  vector,  s  >  0  error  control  threshold,  e  e  (0, 1)  decrease 
factor  for  s,  sample  set  Q  =  (v1,v2,...)'  generated  by  independent  sampling  of  V,  and  n 
an  integer. 

Step  0:  Set  i  =  0  and  z0j  =  0 . 

Step  1:  Compute  (xM,z0M)  =  A^fx^zJ 

Step  2:  Compute  0NiPi(xM,zOM)  and  1Av,a(A+i  >  zo/+i ) 

^  (Xi+1  ’  Z()i-\ )  —  an(I  (Xi+\  ’  Z0i+1  )  —  £ 

Set  N  =  min|[5A; J,lxl04| , 

Ni+l  =  Ni+N, 

pm=nmi  iooo, 

and  replace  s  by  £e . 

else 

Ni+l  =  N, ,  pM  =pr 
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Step  3:  Replace  i  by  i  + 1  and  go  to  Step  1. 

2.  Algorithm  5:  Adaptive  Sample  Size  Algorithm  with  Active-Set 
Strategy 

Algorithm  5  is  a  modification  of  Algorithm  4.  In  this  algorithm  we  only  consider 
a  working  set  of  sample  points,  denoted  by  S,  in  each  iteration.  We  use  the  same 
parameters  as  in  Algorithm  4,  but  additionally  introduce  a  parameter  y  that  we  use  to 
determine  which  sample  points  to  include  in  S .  Algorithm  5  is  identical  to  Algorithm  4 
when  y  =  co  .  Let  A'NpS{x,z0)  denote  the  solution  from  n  iterations,  starting  from  x  and 

z0 ,  of  the  solver  as  applied  to  the  following  problem: 

BP NpS  :  min  f(x ) 

x,z0 

z04 - - - Y^(x,z0)<0,  (49) 

xeX. 

We  also  define  0NpS(x,zQ )  by  (48)  where  (46)  is  replaced  by 

¥Nps (A z0 )  =  max \gNpS {x, z0 ),  max  fJ (x)J ,  (50) 

where 

?NpS  (x’  Zo)  =  Z0  +  —  YsiM.  (51) 

Algorithm  5: 

Data:  As  in  Algorithm  4  and  y  >  0  . 

Step  0:  Set  i  =  0,  z0i  =0,  and  S.  =  {v^v2,...,!^'} . 

Step  1:  Compute  (xM,z()M)  =  A^ix^zJ  and 
S  =  {i  e  {1, 2, ..., Nt}  |  g(xM ,  vJ )  +  y  >  o} . 

Step  2:  Compute  0N^(xi+l,zOM)  and  if/^(xM,z0M) 

If  °xiP/xm’zom)^-£  and  i//Np/xM,z0M)<  £ 
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Set  N- min|[5jV;. J,lxl04| , 

NM=Ni+N, 

SM  ={je{  1,2,...,  Nm}  |  g(xM ,  vJ )  +  y  >  O} 
pM  =  Ni+l  / 1000 ,  and  replace  s  by  se . 

else 

=  Nt ,  pM  =  pt ,  and  SM  =  S,  u  S. 

Step  3:  Replace  i  by  i  + 1  and  go  to  Step  1. 

We  next  compare  the  five  algorithms  described  in  this  chapter  on  six  test 
examples  from  the  literature. 
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IV.  NUMERICAL  EXAMPLES 


In  order  to  evaluate  the  quality  of  solutions  obtained  by  Algorithms  1-5  we  use  a 
method  proposed  by  Royset  (2010a).  That  procedure  gives  a  confidence  interval  for  a 
measure  of  distance,  called  theta,  between  a  design  x  and  a  Fritz-John  point  (e.g., 
Bertsekas,  1999,  pp.  323-335)  of  BP  and  a  confidence  interval  for  constraint  violation  in 
BP.  If  the  confidence  interval  for  theta  for  a  given  x  degenerates  to  the  point  zero  and  the 
constraint  violation  is  nonpositive,  then  x  is  a  feasible  Fritz-John  point  of  BP.  If  the  left 
end  point  of  the  confidence  interval  for  theta  as  well  as  the  right  end  point  of  the 
confidence  interval  for  the  constraint  violation  are  close  to  zero  for  a  specific  x ,  then  x 
is  near  a  Fritz-John  point. 

We  test  Algorithms  1-5  on  six  engineering  design  examples  from  the  litrature. 
These  examples  range  from  simple  models  with  two  design  variables  to  complicated 
structural  designs;  see  Table  2  for  an  overview  of  the  examples  with  number  of  design 
variables  (#  DV),  limit  state  functions  (#  LS),  and  random  variables  (#  RV).  The 
examples  include  a  mix  of  both  linear  and  nonlinear  objective  and  limit-state  functions. 
We  refer  to  the  Appendix  for  detailed  information  about  the  examples.  We  implement 
Algorithms  1-5  in  MATLAB  and  use  the  TOMLAB/SNOPT  (Holmstrom,  1999)  as 
nonlinear  solver.  The  computations  are  run  on  a  desktop  computer  with  3.25  GB  RAM 
and  3.16  GHz  processor  speed.  We  use  1  -  a  =  0.001349898  in  examples,  which 
corresponds  to  the  -3  quantile  of  the  standard  normal  distribution. 
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Example 

Description 

#  DV 

#  LS 

#  R V 

1 

Analytical  example  (Hock  &  Schittkowski, 

1981). 

2 

2 

2 

2 

Design  of  a  cantilever  beam  subject  to  horizontal 

and  vertical  loads  (Eldred  &  Binchon,  2006). 

2 

2 

4 

3 

Design  of  a  rectangular  short  column  (Bichon, 

Mahadevan,  &  Eldred,  2009). 

2 

1 

3 

4 

Design  of  a  uniform  column  of  tubular  section 

with  hinge  joints  at  both  ends  subject  to  a 

random  compressive  load  (Rao,  2009,  pp.  10- 

14). 

2 

2 

1 

5 

Design  of  a  speed  reducer  (Rao,  2009,  pp.  472- 

473). 

7 

9 

7 

6 

Vehicle  design  considering  side-impact 

crashworthiness  (Samson  et  al.  2009). 

11 

10 

7 

Table  2.  Overview  of  Examples  (#  DV  Denotes  Number  of  Design  Variables,  #  LS 
Denotes  Number  of  Limit-state  Functions,  #  RV  Denotes  Number  of  Random 

Variables) 


Tables  3-8  show  numerical  results  for  Algorithms  1-3  when  applied  to  Examples 
1-6.  Here,  we  set  algorithm  parameters  e,  n,  p,  and  v  equal  to  0.001,  5,  1000,  and 

lxlO-6  respectively.  The  sample  size  N  is  varied.  We  stop  Algorithms  1  and  3  when 
SNOPT’s  default  optimality  tolerance  is  satisfied.  Algorithm  2  stops  as  specified  by  its 
Step  2. 
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Example  1 

Algorithm 

Objective 

Function  Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

15.867436 

5.9 

(-0.0901,0] 

N=  1000 

2 

15.867436 

0.1 

(-0.0905,0] 

3 

15.867808 

1.6 

(-0.0779,0] 

1 

15.873157 

213.2 

(-0.0163,0] 

N  =  5000 

2 

15.873157 

0.6 

(-0.0165,0] 

3 

15.873291 

7.4 

(-0.0217,0] 

1 

15.872825 

1406.1 

(-0.0042,0] 

N=  10000 

2 

15.872824 

0.7 

(-0.0032,0] 

3 

15.873006 

15.9 

(-0.0027,0] 

Table  3.  Results  for  Algorithms  1-3  on  Example  1 


Table  3  shows  that  Algorithms  1-3  with  a  larger  sample  size  yield  a  smaller  theta 
interval,  which  means  that  the  corresponding  solutions  become  closer  to  a  Fritz-John 
point  as  N  increases.  That  is,  the  design  quality  improves  as  N  increases.  However,  the 
solution  time  for  Algorithm  1  increases  dramatically  with  larger  sample  sizes. 


Example  2 

Algorithm 

Objective 

Function  Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

9.592329 

10.6 

(-724.3833,0] 

N=  1000 

2 

9.592329 

0.4 

(-723.3969,0] 

3 

9.592329 

24.4 

(-707.6998,0] 

1 

9.697515 

298.7 

(-379.4498,0] 

N  =  5000 

2 

9.697515 

48.9 

(-385.1690,0] 

3 

9.697515 

30.3 

(-384.0893,0] 

1 

9.818805 

1966.6 

(-230.7896,0] 

N=  10000 

2 

9.818805 

199.4 

(-225.4224,0] 

3 

9.818805 

42.3 

(-228.1197,0] 

Table  4.  Results  for  Algorithms  1-3  on  Example  2 


Table  4  shows  that  the  designs  computed  with  the  given  sample  sizes  are  not 
particularly  close  to  a  Fritz-John  point  as  the  theta  intervals  are  large.  For  this  example 
we  compute  a  design  with  theta  interval  (-1.3523,0]  with  a  sample  of  size  2xl05  in 
1207.8  seconds  using  Algorithm  3. 
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Example  3 

Algorithm 

Objective 

Function  Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

236.350524 

5.3 

(-0.0468,0] 

N=  1000 

2 

236.350525 

0.1 

(-0.0475,0] 

3 

236.422249 

4.9 

(-0.0434,0] 

1 

246.780726 

152.2 

(-0.0452,0] 

N  =  5000 

2 

246.780726 

0.2 

(-0.0474,0] 

3 

246.794029 

13.7 

(-0.0463,0] 

1 

247.759968 

1020.8 

(-0.0644,0] 

N=  10000 

2 

247.760098 

1.1 

(-0.0657,0] 

3 

247.769011 

30.3 

(-0.0636,0] 

Table  5. 

Results  for  Algorithms  1-3  on  Example  3 

Example  4 

Algorithm 

Objective 

Function  Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

26.726018 

2.0 

(-0.6408,0] 

N=  1000 

2 

26.726018 

0.2 

(-0.6388,0] 

3 

26.726054 

5.1 

(-0.6323,0] 

1 

26.740881 

17.3 

(-0.0916,0] 

N  =  5000 

2 

26.740881 

2.8 

(-0.0916,0] 

3 

26.740904 

31.6 

(-0.0982,0] 

1 

26.742723 

435.8 

(-0.0982,0] 

N=  10000 

2 

26.742723 

9.8 

(-0.0982,0] 

3 

26.742745 

56.9 

(-0.0982,01 

Table  6. 

Results  for  Algorithms  1-3  on  Example  4 

Example  5 

Algorithm 

Objective 

Function  Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

3469.238339 

3.2 

(-8.2179,0] 

N=  1000 

2 

3469.238338 

0.2 

(-8.2447,0] 

3 

3472.067164 

14.1 

(-8.2455,0] 

1 

3723.733266 

17.8 

(-0.7182,0] 

N  =  5000 

2 

3723.733267 

1.1 

(-0.6977,0] 

3 

3725.084942 

86.9 

(-0.7243,0] 

1 

3760.540797 

1500.0 

(-0.1234,0] 

N=  10000 

2 

3760.540798 

4.3 

(-0.1256,0] 

3 

3761.796161 

156.1 

(-0.1216,0] 

Table  7.  Results  for  Algorithms  1-3  on  Example  5 
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We  see  in  Tables  3-7  that  for  Algorithm  2,  the  solution  times  with  different 
sample  sizes  are  close,  meaning  that  Algorithm  2  is  not  highly  affected  by  the  increase  in 
sample  size. 


Example  6 

Algorithm 

Objective 

Function  Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

24.597346 

2.6 

(-0.4643,0] 

N=  1000 

2 

24.597347 

0.2 

(-0.4619,0] 

3 

24.605073 

7.3 

(-0.4532,0] 

1 

24.774812 

20.8 

(-0.0288,0] 

N  =  5000 

2 

24.774808 

0.8 

(-0.0289,0] 

3 

24.779162 

29.6 

(-0.0265,0] 

1 

24.850419 

118.0 

(-0.0124,0] 

N=  10000 

2 

24.850436 

1.9 

(-0.0135,0] 

3 

24.853971 

68.8 

(-0.0086,01 

Table  8.  Results  for  Algorithms  1-3  on  Example  6 


Tables  3-8  show  that  the  objective  function  values  computed  with  Algorithms  1 
and  2  are  close  to  each  other.  However,  the  solution  times,  which  are  quite  similar  for  small 
sample  sizes,  differ  as  we  increase  the  sample  size.  This  is  because  we  add  one  more 
variable  and  constraint  for  every  limit-state  function  for  each  increment  in  the  sample  size. 
For  example,  the  Example  5  solution-time  ratio  for  Algorithm  1  to  Algorithm  2  is  16  for 
N  — 1000 ,  while  it  is  349  for  N  =  10000.  Therefore,  we  conclude  that  for  more  complex 
problems,  Algorithm  1  has  large  memory  and  longer  solution  time  requirements.  For 
Algorithm  3,  the  objective  function  value  is  worse  than  that  of  Algorithms  1  and  2. 
However,  we  do  not  see  a  dramatic  increase  in  solution  times  with  larger  sample  sizes  for 
Algorithm  3.  We  have  the  largest  increase  in  solution  time  in  Example  4,  where  the  time 
for  N  =  10000  is  1 1  times  larger  than  that  for  N  =  1000.  However,  the  solution  time  for 
Algorithm  1  increases  by  a  factor  of  218  for  the  same  example.  We  next  discuss  the 
differences  in  solutions  for  different  algorithm  parameter  changes. 

Although  we  see  that  Algorithm  2  is  able  to  solve  Examples  1-6  quickly  and  with 
high  accuracy,  the  choice  of  user-defined  algorithm  parameters  yields  differences  in 
solution  times.  Tables  9-12  present  solution  times  for  Algorithm  2  in  seconds  for  different 
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sample  sizes  as  the  algorithm  parameters  change.  The  row  numbers  0.01,  0.001  and  0.0001 
and  the  column  numbers  from  1  to  10  are  the  values  for  s  and  n ,  respectively. 
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Iteration  Limit  for  the  Solver  in) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0.01 

28.7 

16.8 

8.4 

5.7 

4.4 

4.8 

4.4 

5.6 

6.2 

8.4 

0.001 

29.2 

17.3 

9.6 

7.9 

6.3 

6.7 

7.7 

7.1 

7.6 

9.6 

0.0001 

29.9 

18.2 

10.9 

8.7 

7.0 

7.8 

7.5 

7.8 

8.3 

10.0 

Table  9.  Run  Times  (sec)  for  Algorithm  2  on  Example  5  with  N  =  12000  Given  Different 

Algorithm  Parameter  Settings 


Table  9  shows  that  as  e  increases  the  solver  needs  more  time.  Moreover,  the 
computation  times  improve  from  n  —  1  through  n  —  5 ,  and  gradually  deteriorate  for 
larger  values  of  n .  But  in  Table  9  we  see  that  the  choice  of  algorithm  parameters  does 
not  yield  much  different  solution  times  for  Example  5  with  N  -  12000 . 
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Iteration  Limit  for  the  Solver  in 

) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0.01 

335 

212 

146 

153 

187 

215 

148 

662 

1157 

149 

0.001 

372 

238 

357 

335 

429 

632 

691 

775 

1139 

1586 

0.0001 

417 

302 

632 

538 

612 

769 

862 

944 

1207 

1691 

Table  10.  Run  Times  (sec)  for  Algorithm  2  on  Example  5  with  N  =  120000  Given 

Different  Algorithm  Parameter  Settings 


We  see  in  Table  10  that  with  a  larger  sample  the  variability  in  solution  times  is 
more  evident.  Moreover,  we  see  that  the  results  in  Table  10  do  not  follow  the 
characteristics  of  those  on  Table  9,  where  the  run  times  are  unevenly  increasing  or 
decreasing  with  respect  to  parameter  choices. 
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Iteration  Limit  for  the  Solver  (n 

) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0.01 

6.24 

2.22 

2.61 

3.44 

3.06 

3.37 

3.08 

3.03 

3.05 

2.99 

0.001 

4.17 

4.65 

5.09 

5.08 

5.12 

5.21 

5.21 

5.28 

5.23 

5.36 

0.0001 

6.43 

5.18 

6.14 

6.11 

6.09 

6.20 

6.19 

6.24 

6.22 

6.19 

Table  1 1.  Run  Times  (sec)  for  Algorithm  2  on  Example  6  with  N  =  12000  Given  Different 

Algorithm  Parameter  Settings 


28 


We  see  in  Table  1 1  that  larger  s  values  yield  longer  run  times.  However,  we  do 
not  see  a  clear  increasing  or  decreasing  pattern  in  solution  times  depending  on  n. 


s 

Iteration  Fimit  for  the  Solver  (n 

) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0.01 

174 

85 

45 

45 

46 

47 

47 

47 

49 

49 

0.001 

117 

130 

191 

150 

152 

151 

149 

151 

150 

150 

0.0001 

254 

199 

305 

284 

288 

287 

287 

287 

287 

288 

Table  12.  Run  Times  (sec)  for  Algorithm  2  on  Example  6  with  N  -  120000  Given 

Different  Algorithm  Parameter  Settings 


Tables  9-12  show  that  there  is  no  single  parameter  value  that  is  best  in  all 
examples  for  Algorithm  2.  Thus,  we  conclude  that  the  parameter  choice  gives  different 
solution-time  results  not  only  for  the  different  examples  but  also  for  the  same  example 
with  different  sample  sizes.  Furthermore,  s  values  larger  than  0.01  effectively  lead  to 
large  working  sets  and  turn  Algorithm  2  into  Algorithm  1 .  We  next  discuss  the  effects  of 
parameter  choice  for  Algorithm  2  on  solution  quality. 

We  select  s  =  0.01  and  n  =  1  as  the  base  parameter  choices  for  Tables  13  and  14. 
We  then  optimize  with  different  values  of  s  and  n,  and  then  calculate  the  absolute 
differences  in  the  objective  function  values  from  that  obtained  using  the  base  values  of  s 
and  n. 


S 

Iteration  Fimit  for  the  Solver  (n 

) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0.01 

0 

0.8 

0.8 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.001 

136 

0.04 

25 

3360 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.0001 

113 

4.2 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

Table  13.  Objective  Function  Value  Differences  (in  1 0‘7)  for  the  Parameter  Selections  with 
Respect  to  Objective  Function  of  the  Base  Parameter  Choice  (with  Algorithm  2 

on  Example  5  with  N  - 12000) 
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£ 

Iteration  Limit  for  the  Solver  (n 

) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

0.01 

0 

3.7 

4.1 

3.7 

3.7 

3.7 

3.7 

3.7 

3.7 

3.7 

0.001 

3.0 

6.7 

5.3 

7.7 

3.7 

3.7 

3.7 

3.7 

3.7 

3.7 

0.0001 

16 

204 

11 

11 

3.7 

3.7 

3.7 

3.7 

3.7 

3.7 

Table  14.  Objective  Function  Value  Differences  (in  1 0‘7)  for  the  Parameter  Selections  with 
Respect  to  Objective  Function  of  the  Base  Parameter  Choice  (with  Algorithm  2 

on  Example  6  with  N  - 120000) 


We  see  in  Tables  13  and  14  that  there  is  no  dramatic  variability  in  objective  function 
values  computed  with  various  parameters.  Thus,  we  conclude  that  the  parameter  choice 
does  not  affect  the  quality  of  the  design.  We  next  discuss  Algorithm  3  with  various 
smoothing  parameter  p  values. 

We  apply  Algorithm  3  on  Example  1  with  a  sample  size  of  15000.  We  optimize 
the  design  with  different  values  of  p,  compute  the  theta  interval  in  order  to  evaluate  the 
design  quality,  and  display  the  results  in  Table  15. 


P 

Objective  Function 

Value 

Solution  Time 

(seconds) 

Theta  Interval 

1 

20.62291723 

24.5 

(-1.4575161,0] 

10 

16.29547551 

17.2 

(-0.5836538,0] 

100 

15.88688647 

27.2 

(-0.3088083,0] 

1000 

15.87386166 

30.8 

(-0.0062922,0] 

10000 

15.87370318 

33.8 

(-0.0064495,0] 

100000 

15.87370048 

35.4 

(-0.0062150,0] 

1000000 

15.87370032 

38.5 

(-0.0063207,0] 

10000000 

15.87370031 

41.1 

(-0.0061170,01 

Table  15.  Algorithm  3  on  Example  1  with  Different  p  Values 


Table  15  shows  that  design  quality  improves  as  we  increase  the  value  of  p . 
However,  we  do  not  have  much  improvement  in  the  theta  interval  after  p  =  1 000 . 
Instead,  we  see  only  run-time  increases  with  larger  values  of  p .  Since  large  p  values 
may  result  in  ill-conditioning  and  slow  convergence  of  the  solver,  we  use  p  =  1000  as 
default. 
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We  conclude  that  Algorithm  1  cannot  handle  large  sample  sizes  due  to  memory 
difficulty.  In  fact,  Algorithm  1  solves  BPV  that  include  more  than  Nm  nonlinear 

constraints  and  N  variables.  Algorithm  2  may  need  to  consider  the  same  number  of 
constraints  and  variables,  but  the  numerical  tests  indicate  that  its  active-set  strategy  is 
reasonably  efficient  in  practice.  On  solutions  of  Examples  1-6,  we  find  that  BP N(W) 

contains  on  average  0.2%  of  nonlinear  constraints  that  we  have  in  BPV .  Finally,  in 

Algorithm  3  in  contrast  to  Algorithms  1  and  2  we  have  only  one  nonlinear  constraint 
associated  with  limit-state  functions  and  number  of  design  variables  does  not  increase 
with  sample  size.  Therefore,  Algorithms  2  and  3  efficiently  compute  near-optimal 
designs  for  large  N .  Algorithm  2  can  handle  large  samples,  but  lacks  the  ability  to  deal 
with  extremely  large  problems.  But  we  apply  Algorithm  3  on  Example  6  with  a  sample 
size  of  lxlO7  and  compute  the  design,  with  a  (-0.0006,0]  theta  interval  in  18  hours. 
Thus,  we  conclude  that  Algorithm  3  is  able  to  solve  BPat  with  exceptionally  large  sample 
sizes.  We  next  discuss  consequences  of  various  algorithm  parameters  for  Algorithms  4 
and  5. 

We  apply  Algorithm  4  to  Example  6  to  see  the  theta  interval  differences  due  to 
the  choice  of  algorithm  parameters.  Table  16  presents  the  results  for  Algorithm  4  with  a 
run-time  limit  used  as  the  stopping  criterion.  That  is,  the  algorithm  stops  after  a  specific 
time,  here  1 5  minutes,  and  evaluates  the  quality  of  the  last  computed  design. 


e 

5 

n 

Theta  Interval 

Final  Sample  Size 

0.1 

0.3 

15 

(-0.0273379,0] 

10598 

0.1 

0.3 

20 

(-0.0175738,0] 

13777 

0.1 

0.5 

15 

(-0.0089000,0] 

45624 

0.1 

0.5 

20 

(-0.0089010,0] 

45624 

0.3 

0.3 

15 

(-0.0202027,0] 

30267 

0.3 

0.3 

20 

(-0.0079490,0] 

69347 

0.3 

0.5 

15 

(-0.0599124,0] 

45624 

0.3 

0.5 

20 

(-0.0041565,0] 

85624 

0.5 

0.3 

15 

(-0.0116882,0] 

89342 

0.5 

0.3 

20 

(-0.0035743,0] 

69347 

0.5 

0.5 

15 

(-0.0054132,0] 

95624 

0.5 

0.5 

20 

(-0.0077265,0] 

95624 

Table  16.  Parameter  Comparison  for  Example  6  with  Algorithm  4  Starting  From  Initial 

N  —  1000  With  15 -Minute  Run  Time  Limit 
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From  Table  16,  we  see  that  for  larger  values  of  e  the  sample  size  tends  to  increase 
more  compared  to  smaller  e  values,  which  in  turn  may  cause  memory  difficulty  for  more 
complex  structural  design  examples.  We  obtain  better  results  with  n  =  20  than  with 
n  -  15  for  fixed  e  and  s  .  We  also  see  that  we  get  generally  better  results  with  s  -  0.5 
when  the  other  parameters  are  hold  constant.  Therefore,  we  conclude  that  e  =  0. 1 , 
s  =  0.5  ,  and  n  =  20  is  an  appropriate  choice  for  parameter  values  and  use  them  in  the 
following  calculations.  We  next  present  Algorithm  5  results  for  different  y  values. 

As  discussed  in  Chapter  III,  we  use  y  to  determine  the  active  sample  points.  We 
apply  Algorithm  5  to  Example  1  with  various  /-values  and  construct  Table  17  after  15- 
minute  runs. 


r 

Constraint  Violation 

Theta  Interval 

Final  Sample  Size 

0.0000001 

(-00,1.52274) 

(-  1.3017056,0] 

75624 

0.000001 

(-00,-0.57666) 

(-3.0074925,0] 

55624 

0.00001 

(-oo  ,0.00037) 

(-0.0041879,0] 

65624 

0.0001 

( -oo  ,0.00009) 

(-0.0017403,0] 

55624 

0.001 

( -co  ,-0.00201) 

(-0.0026550,0] 

45624 

0.01 

(-oo  ,0.00037) 

(-0.0004014,0] 

45624 

1 

( -co  ,-0.00028) 

(-0.0001733,0] 

145624 

10 

( -co  ,-0.00024) 

(-0.0009416,0] 

385624 

100 

( -oo  ,0.00047) 

(-0.0011301,0] 

425624 

1000 

( -co  ,-0.00 123) 

(-0.0009869,0] 

95624 

10000 

( -oo  ,0.00056) 

(-0.0007012,0] 

95624 

100000 

( -oo  ,0.00050) 

(-0.0009243,0] 

105624 

1000000 

(-oo  ,-0.000 14) 

(-0.0019971,0] 

125624 

10000000 

( -oo  ,0.00026) 

(-0.0018027,0] 

135624 

Table  17.  y  Comparison  With  Algorithm  5  on  Example  1  with  Initial 

A  =  1000,  e  =  0.1,  5  =  0.5,  n  =  20  with  15-Minute  Run  Time  Limit 


Table  17  shows  that  solution  quality  improves  as  we  increase  y  from  0.0000001 
to  1,  but  degrades  for  y  >  1 .  Since  the  smallest  theta  interval  appears  to  occur  at  y  =  l, 
we  want  to  investigate  the  values  of  y  near  1,  and  use  0.1,  1,  and  10,  in  the  following 
calculations. 

We  present  below  the  results  for  Algorithms  4  and  5  on  Examples  1-6.  In  order  to 
show  the  performance  of  the  algorithms  we  set  a  run-time  limit  as  the  stopping  criterion. 

We  use  different  run  times  in  accordance  with  the  complexity  of  the  example.  For 
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relatively  simple  Examples  1  and  3,  we  set  the  time  limit  to  15  and  30  minutes, 
respectively;  the  time  limit  is  set  to  60  minutes  for  the  rest.  Rather  than  setting  a  time 
limit,  however,  the  user  may  prefer  to  modify  Algorithms  4  and  5  to  iteratively  check  the 
design  quality  and  stop  when  a  user-defined  quality  level  is  reached.  We  start  with  initial 
N  = 1000 . 


e 

5 

n 

Y 

Final  Sample 

Size 

Constraint 

Violation 

Theta  Interval 

Alg.4 

0.1 

0.5 

5 

- 

45624 

(-00  ,0.000 14) 

(-0.00036,0] 

0.3 

0.5 

5 

- 

105624 

( -co, -0.00 104) 

(-0.00047,0] 

0.1 

0.5 

20 

0.1 

145624 

(-oo, -0.00028) 

(-0.00017,0] 

Alg.5 

0.1 

0.5 

20 

1 

385624 

(-oo, -0.00024) 

(-0.00094,0] 

0.1 

0.5 

20 

10 

425624 

( -oo  ,0.00046) 

(-0.00113,01 

Table  18.  Results  for  Algorithms  4  and  5  on  Example  1  with  15-Minute  Run  Time 


We  see  in  Table  18  that  we  have  the  best  theta  interval  with  Algorithm  5,  when 
y  =  0.1 .  The  design  computed  by  Algorithm  4  with  e  —  0.1  has  also  a  small  theta  interval 
and  its  constraint  violation  is  less  than  for  Algorithm  5 . 


e 

5 

n 

Y 

Final  Sample 

Size 

Constraint 

Violation 

Theta  Interval 

Alg.4 

0.1 

0.5 

20 

- 

335624 

(-oo,0.94085) 

(-1.24864,0] 

0.3 

0.5 

20 

- 

205624 

(-oo,1.79601) 

(-2.15182,0] 

0.1 

0.5 

20 

0.1 

194581 

(-oo,0.09473) 

(-0.54203,0] 

Alg.5 

0.1 

0.5 

20 

1 

2891871 

(-oo,0.19126) 

(-0.66174,0] 

0.1 

0.5 

20 

10 

3491871 

(-oo,0.40577) 

(-0.60458,0] 

Table  19.  Results  for  Algorithms  4  and  5  on  Example  2  with  60-Minute  Run  Time 


Table  19  shows  that  the  sample  size  rapidly  increases  for  Algorithm  5.  However, 
Algorithm  5  only  considers  active  constraints.  Therefore,  Algorithm  4  is  more  likely  to 
cause  memory  difficulty  for  highly  complex  structural  design  examples. 
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e 

5 

n 

Y 

Final  Sample 

Size 

Constraint 

Violation 

Theta  Interval 

Alg.4 

0.1 

0.5 

20 

- 

215624 

( -oo  ,0.00000) 

(-0.00863,0] 

0.3 

0.5 

20 

- 

205624 

(-oo,0.00123) 

(-0.00676,0] 

0.1 

0.5 

20 

0.1 

691871 

(-oo,7.43356) 

(-1.67174,0] 

Alg.5 

0.1 

0.5 

20 

1 

1591871 

( -oo  ,0.00034) 

(-0.00941,0] 

0.1 

0.5 

20 

10 

2491871 

( -oo  ,0.00095) 

(-0.00437,01 

Table  20.  Results  for  Algorithms  4  and  5  on  Example  3  with  30-Minute  Run  Time 


From  Table  20,  we  see  that  Algorithm  5  with  y  =  0.1  results  in  an  unsatisfactory 
solution.  Since  the  theta  interval  is  large  and  the  right  end  point  of  the  constraint  violation 
interval  is  also  large,  we  conclude  that  the  solution  is  infeasible. 


e 

5 

n 

Y 

Final  Sample 

Size 

Constraint 

Violation 

Theta  Interval 

Alg.4 

0.1 

0.5 

20 

- 

85624 

(-oo,-0.04217) 

(-0.09810,0] 

0.3 

0.5 

20 

- 

75624 

(-co,0.07769) 

(-0.18332,0] 

0.1 

0.5 

20 

0.1 

119721 

( -oo  ,0.00000) 

(-0.17360,0] 

Alg.5 

0.1 

0.5 

20 

1 

129721 

(-co,0.52574) 

(-0.65265,0] 

0.1 

0.5 

20 

10 

991871 

(-oo,-0.02971) 

(-0.09809,0] 

Table  2 1 .  Results  for  Algorithms  4  and  5  on  Example  4  with  60-Minute  Run  Time 


Table  21  shows  that  for  Example  4,  the  quality  of  the  designs  by  Algorithm  4  with 
e  =  0. 1  and  Algorithm  5  with  y  =  10  are  almost  equal. 


e 

5 

n 

Y 

Final  Sample 

Size 

Constraint 

Violation 

Theta  Interval 

Alg.4 

0.3 

0.5 

20 

- 

11389 

(-oo,1.85124) 

(-1.96464,0] 

0.8 

0.5 

20 

- 

65624 

(-oo,0.21324) 

(-0.24076,0] 

0.8 

0.5 

20 

0.1 

129721 

(-oo,0.21682) 

(-0.23135,0] 

Alg.5 

0.8 

0.5 

20 

1 

129721 

(-oo  ,0.04015) 

(-0.05877,0] 

0.8 

0.5 

20 

10 

129721 

(-oo,0.08343) 

(-0.10257,01 

Table  22.  Results  for  Algorithms  4  and  5  on  Example  5  with  60-Minute  Run  Time 


We  see  in  Table  22  that  for  Example  5,  Algorithm  5  yields  better  theta  intervals 
and  small  constraint  violations  than  Algorithm  4. 


34 


e 

5 

n 

Y 

Final  Sample 

Size 

Constraint 

Violation 

Theta  Interval 

Alg.4 

0.1 

0.5 

20 

- 

65624 

(-°°  ,0.00659) 

(-0.01540,0] 

0.3 

0.5 

20 

- 

125624 

(-<»  ,0.00000) 

(-0.00464,0] 

0.1 

0.5 

20 

0.1 

2491871 

(-°°  ,0.00003) 

(-0.00168,0] 

Alg.5 

0.1 

0.5 

20 

1 

2691871 

(-°°  ,0.00041) 

(-0.00183,0] 

0.1 

0.5 

20 

10 

2591871 

(-00  ,0.00029) 

(-0.00068,0] 

Table  23.  Results  for  Algorithms  4  and  5  on  Example  6  with  60-Minute  Run  Time 


From  Table  23,  we  see  that  we  obtain  a  good  design  with  Algorithm  5  when 
/  =  10.  The  lower  end  point  of  the  theta  interval  is  almost  the  same  as  that  obtained  after 

1 8  hours  using  Algorithm  3  with  N  =  1  x  1 07 . 

We  see  in  Tables  18-23  that  Algorithm  5  is  capable  of  handling  larger  sample 
sizes,  and  thus  Algorithm  5  generally  computes  better  designs  compared  to  Algorithm  4. 
Moreover,  because  Algorithm  5  only  considers  the  active  sample  points  the  iterations 
take  less  time  compared  to  Algorithm  4.  Therefore,  with  a  stopping  criterion  based  on  a 
user-defined  quality  level,  Algorithm  5  is  likely  to  be  faster. 

We  next  present  failure  probability  and  buffered  failure  probability  comparisons 
for  a  given  design.  For  each  example,  we  select  the  design  with  the  closest  theta  interval 
among  the  solutions  presented  in  Tables  3-8  and  Tables  18-23.  We  calculate  a  95% 
confidence  interval  (Cl)  using  MCS  for  the  failure  probability.  Using  the  same  sample, 
we  also  estimate  the  buffered  failure  probability.  Since  we  use  the  same  sample,  we 
assume  that  the  error  in  the  buffered  failure  probability  calculation  is  similar  to  that  in  the 
failure  probability  estimate.  Table  24  presents  the  resulting  estimates.  We  see  that  the 
buffered  failure  probability  is  greater  than  the  failure  probability  for  the  same  design  as 
expected;  see  Chapter  I.  On  average,  the  buffered  failure  probability  is  30%  of  the  failure 
probability. 
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Example 

Design 

Theta 

Interval 

Cl  for  Failure  Probability 

Buffered  Failure 

Probability 

Estimate 

1 

8.90895 

2.81726 

(-0.00017,0] 

(0.00052,+/-  0.00005) 

0.001329 

2 

3.88551 

2.73362 

(-0.54202,0] 

(0.000001,+/-  0.0000006) 

0.001015 

3 

9.82582 

25 

(-0.00437,0] 

(0.00052,+/-  0.00005) 

0.001402 

4 

5.45094 

0.29593 

(-0.09820,0] 

(0.00037,+/- 0.000036) 

0.000971 

5 

3.6 

0.72 

19.52866 

7.56277 

8.28022 

3.47997 

5.40634 

(-0.12156,0] 

(0.00047, +/- 0.000046) 

0.004488 

6 

0.5 

1.43498 

0.5 

1.25441 

1.05796 

0.5 

0.34 

0.345 

15 

15 

(-0.00068,0] 

(0.00031,+/- 0.00003) 

0.001382 

Table  24.  Computational  Comparisons  of  Failure  Probabilities  and  Buffered  Failure 

Probabilities 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

Reliability-based  design  optimization  (RBDO)  aims  to  detennine  a  minimum-cost 
design  for  an  engineering  structure  subject  to  one  or  more  reliability  constraints;  this 
thesis  considers  models  having  a  single  reliability  constraint.  Traditionally,  the  reliability 
constraint  is  given  in  terms  of  an  upper  bound  on  the  failure  probability,  i.e.,  the 
probability  that  the  system  performs  unsatisfactory.  This  approach  is  theoretically  and 
computationally  troublesome  since  a  constraint  on  the  failure  probability  is  difficult  to 
deal  with  in  the  algorithms  used  to  solve  the  nonlinear  programs  arising  in  RBDO.  This 
thesis  considers  an  alternative  approach  to  RBDO  based  on  “buffered  failure  probability,” 
and  examines  five  solution  algorithms,  four  of  which  are  developed  in  this  thesis,  that  use 
sample-average  approximations.  Buffered  failure  probability  is  more  conservative  than 
the  traditional  failure  probability.  Thus,  a  design  that  satisfies  a  reliability  constraint 
based  on  the  buffered  failure  probability  also  satisfies  one  based  on  the  failure 
probability.  Finally,  the  buffered  failure  probability  is  much  easier  to  handle  in 
optimization  algorithms  as  it  results  in  smooth  and  possibly  convex  optimization 
problems. 

The  buffered  failure  probability  approach  uses  sample  averages  to  estimate 
expectations.  In  numerical  tests,  we  show  that  the  sample  size  needs  to  be  set  relatively 
large  to  ensure  high-quality  solutions,  and  that  one  standard  algorithm  may  break  down 
because  of  the  resulting  memory  and  run-time  requirements.  We  develop  new  algorithms 
that  overcome  this  difficulty  and  obtain  an  average  speed-up  in  solution  time  by  a  factor 
of  560  in  comparison  with  the  existing  methodology  based  on  a  standard  nonlinear 
solver.  We  are  able  to  handle  sample  sizes  two  orders  of  magnitude  larger  in  comparison 
with  the  existing  method.  We  also  avoid  the  need  for  preselecting  sample  sizes,  which 
can  be  difficulty  in  practice,  by  using  adaptive  sample-adjustment  schemes. 

We  examine  the  difference  between  the  failure  probability  and  buffered  failure 
probability  approaches  in  a  stochastic  knapsack  problem  as  well  as  other  examples  and 
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find  that  for  these  test  examples  the  buffered  failure  probability  averages  typically  three 
times  larger  than  the  failure  probability  for  a  design  computed  with  buffered  failure 
probability  approach.  Hence,  a  design  based  on  the  buffered  failure  probability  approach 
may  result  in  a  more  costly  but  more  reliable  design  than  that  obtained  using  a  failure 
probability  approach.  However,  in  view  of  the  computational  difficulties  associated  with 
this  approach,  we  believe  that  the  buffered  failure  probability  approach  is  a  viable 
alternative. 

B.  SUGGESTED  WORK  AHEAD 

This  research  area  is  open  to  more  developments  in  efficiency  and  accuracy  of  the 
algorithms.  We  believe  the  active-set  strategy  deserves  the  interest  of  future  researchers 
with  the  goal  to  eliminate  its  sensitivity  to  user-specified  parameters. 
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APPENDIX 


This  appendix  describes  the  six  computational  examples  tested  in  this  thesis. 


Example  1:  Analytical  problem  (Hock  &  Schittkowski,  1981). 


DESIGN  VARIABLES 

Variable 

Description 

Lower 

bound 

Upper 

bound 

A 

Design  variable 

2 

50 

*2 

Design  variable 

0 

50 

Table  25.  Design  Variable  Descriptions  and  Bounds  for  Example  1 


RANDOM  VARIABLES 

V 

Description 

Distribution 

Parameters 

(A  v) 

VI 

Random  vector 

Normal 

(25,  0.03) 

V2 

Random  vector 

Normal 

(25,  0.03) 

Table  26.  Random  Variable  Descriptions  and  Distribution  Parameters  for  Example  1 


Objective  function: 
min  0.  IaVV 
Limit-state  functions: 

gx(x,v)  =  vl-xlx2 

g2(x,v)  =  v2-x2]  -x\ 
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Example  2:  Cantilever  beam  (Eldred  &  Binchon,  2006). 


DESIGN  VARIABLES 

Variable 

Description 

Lower 

bound 

Upper 

bound 

Thickness  of  the  beam  (cm) 

1 

4 

x2 

Width  of  the  beam  (cm) 

1 

4 

Table  27.  Design  Variable  Descriptions  and  Bounds  for  Example  2 


RANDOM  VARIABLES 

V 

Description 

Distribution 

Parameters 

(/4<T) 

VI 

Yield  stress 

Normal 

(4xl0\2xl03) 

V2 

Young’s  Modulus 

Normal 

(2.9x107,1.45x106) 

V3 

Horizontal  loads  (kg) 

Lognormal 

(5,0.5) 

V4 

Vertical  loads  (kg) 

Lognormal 

(5,0.5) 

Table  28.  Random  Variable  Descriptions  and  Distribution  Parameters  for  Example  2 

Objective  function: 

min  xxx2 

Limit-state  functions: 


£i(w) 


600 v4  600 v3 


,  4000000 

g,(x,v)  =  - 

v2Xi*2 


2.25 
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Example  3:  Rectangular  short  column  (Bichon,  Mahadevan,  &  Eldred,  2009). 


DESIGN  VARIABLES 

Variable 

Description 

Lower 

bound 

Upper 

bound 

A 

Column  cross  section  width  (cm) 

5 

15 

x2 

Column  cross  section  depth  (cm) 

15 

25 

Table  29.  Design  Variable  Descriptions  and  Bounds  for  Example  3 


RANDOM  VARIABLES 

V 

Description 

Distribution 

Parameters 

(M,  o) 

VI 

Axial  force  (kg) 

Normal 

(500,100) 

V2 

Bending  moment 

Normal 

(2000,400) 

V3 

Yield  stress 

Lognormal 

(5,0.5) 

Table  30.  Random  Variable  Descriptions  and  Distribution  Parameters  for  Example  3 


Objective  function: 

min  xtx2 

Limit-state  function: 


g(x,v) 


4v2 

AUG 


+  - 


2  2  2 
A  UU 


1 
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Example  4:  Tubular  column  (Rao,  2009,  pp.  10-14). 

DESIGN  VARIABLES 


Variable 

Description 

Lower 

bound 

Upper 

bound 

A 

Diameter  of  the  column  (cm) 

2 

14 

x2 

Thickness  of  the  tube  (cm) 

0.2 

0.8 

Table  3 1 .  Design  Variable  Descriptions  and  Bounds  for  Example  4 


RANDOM  VARIABLES 

V 

Description 

Distribution 

Parameters 

0 u,cx ) 

V 

Load  on  the  system  (kg) 

Normal 

(2500,  10) 

Table  32.  Random  Variable  Descriptions  and  Distribution  Parameters  for  Example  4 


Objective  function: 
min  9.82.r1.r2  +  2x] 

Limit-state  functions: 

&(*.v)  =  — - 500 

nxxx2 

g2(x,v)  =  — - \.ln2{x\-x22) 

nx  jX2 
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Example  5:  Speed  reducer  (Rao,  2009,  pp.  472-473). 


DESIGN  VARIABLES 

Variable 

Description 

Lower 

bound 

Upper 

bound 

A 

Face  width  (cm) 

2.6 

3.6 

x2 

Module  of  teeth  (cm) 

0.7 

0.8 

x3 

Number  of  teeth  on  pinion 

17 

28 

x4 

Length  of  shaft  1  (cm) 

7.3 

8.3 

x5 

Length  of  shaft  2  (cm) 

7.3 

8.3 

*6 

Diameter  of  shaft  1  (cm) 

2.9 

3.9 

x7 

Diameter  of  shaft  2  (cm) 

5.0 

5.5 

Table  33.  Design  Variable  Descriptions  and  Bounds  for  Example  5 


RANDOM  VARIABLES 

V 

Description 

Distribution 

Parameters 

(A  c) 

VI 

Material  property 

Normal 

(xr  0.03) 

Material  property 

Normal 

(x2,  0.03) 

Material  property 

Normal 

(jc3,  0.03) 

K 

Material  property 

Normal 

(x4,  0.03) 

Material  property 

Normal 

(x5,  0.03) 

K 

Material  property 

Normal 

(x6,  0.03) 

v7 

Material  property 

Normal 

(x7,  0.03) 
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Table  34.  Random  Variable  Descriptions  and  Distribution  Parameters  for  Example  5 


Objective  function: 

min  0.7854xjX2 (3. 33333x2  +  14.9334x3  -43.0934)-1.508X[(Xg  +  x2)  + 
7.477(Xg  +  x7)  +  0.7854(x4Xg  +  x5x2) 

Limit-state  functions: 


gt(x,v)  = 


27 


ViV22v3 


-1 


S'2(w)  = 


397.5 

ViV22V32 


-1 


g3(w) 


1-93v43 

v2v3v64 


-1 


g4(w)  = 


l-93vg 

v2v3v74 


-1 


g5(w)  = 


A745v4V 
V  Vv3  y 


+  1.69x10 


0.5 


O.lv 


1100 


g6(x’v)  = 


/745v5V 

V  V2V3  y 


+  1. ,575x10s 


0.5 


O.lv, 


850 


g7(x,v)  =  v2v3-40 


g8(w) 


(1-5v6+1-9) 

V4 


g9(w) 


(1 .  lv7  + 1 .9) 
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Example  6:  Vehicle  side-impact  crashworthiness  problem  (Samson  et  al.  2009). 


DESIGN  VARIABLES 

Variable 

Description 

Lower 

bound 

Upper 

bound 

X l 

Thickness  of  B-pillar  inner  (cm) 

0.5 

1.5 

*2 

Thickness  of  B-pillar  reinforce  (cm) 

0.5 

1.5 

x3 

Thickness  of  floor  side  inner  (cm) 

0.5 

1.5 

x4 

Thickness  of  cross  member  (cm) 

0.5 

1.5 

x5 

Thickness  of  door  beam  (cm) 

0.5 

1.5 

*6 

Thickness  of  door  belt  line  (cm) 

0.5 

1.5 

X-J 

Thickness  of  roof  rail  (cm) 

0.5 

1.5 

*8 

Material  property  of  B-pillar  inner 

0.345 

0.345 

Xg 

Material  property  of  floor  side  inner 

0.345 

0.345 

Vo 

Barrier  hitting  height  (cm) 

15 

15 

Vi 

Barrier  hitting  position  (cm) 

15 

15 

Table  35.  Design  Variable  Descriptions  and  Bounds  for  Example  6 


49 


RANDOM  VARIABLES 

V 

Description 

Distribution 

Parameters 

(A  &) 

v; 

Material  property 

Normal 

(xp  0.03) 

v2 

Material  property 

Normal 

(x2,  0.03) 

Material  property 

Normal 

(x3,  0.03) 

Material  property 

Normal 

(x4,  0.03) 

Material  property 

Normal 

(x5,  0.03) 

K 

Material  property 

Normal 

(x6,  0.03) 

v7 

Material  property 

Normal 

(x7,  0.03) 

Table  36.  Random  Variable  Descriptions  and  Distribution  Parameters  for  Example  6 


Objective  functions: 

min  1.98  +  4.9x1  +  6.67 x2  +  6.98x3  +  4.01x4  +  1.78x5  +  2.73xv 
Limit-state  function: 

gi(x,v)  =  1.16  -  0.3717v2v4  -  0.0093  lv2x10  -  0.484v3x9  +  0.01343v6x10  -  1 


g2(x,v)  =  0.261- 0.01  59vjV2  -0.188vjX8  -0.019v2v7  +  0.0144v3v5  +  0.0008757v5x1g  + 
0.080445v6x9  +0.00139xgxn  - 0.00001 575x10xn  -0.32 

g3(x,v)  =  0.2 14  + 0.008  17v5  -0.131vjX8  -0.0704vjX9  +0.03099v2v6  -0.018v2v7  + 

0.0208v3xg  +0.121v3x9  -0.00364v5v6  +  0.00077  15v5x10  -0.0005354v6x10  + 
0.00121x8xn  -0.32 


g4(x,v)  =  0.74-0.61v2  -0.163v3x8  +0.001232v3x10  -0.166v7x9  +0.0227v2  -0.32 


g5(x,v)  =  28.98  -  3.8 lv3  -  4.2v,v2  +  0.0207v5x1G  +  6.63v6x9  -  7.7v7x8  +  0.32x9x10  -  32 
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g6(x,v)  =  33.86  +  2.95v3  +  0.1792x10  -5.057v,v2  -1  lv2x8  -0.0215v5x10  -9.98v7x8  + 

22x8x9  -32 

g7(x,v)  =  46.36  -  9.9v2  -  12.9v,x8  +  0.1 107v3x10  -  32 

g8(x,  v)  =  4.72  -  0.5v4  -  0.19v2v3  -  0.0122v4x10  +  0.009325v6x10  +  0.000191x7,  -  4 
gg(x,v)  =  10.58  -  0.674v,v2  -  1.95v2x8  +  0.02054v3x10  -  0.0198v4x10  +  0.028v6x10  -  9.9 
g10(x, v)  =  16.45  -  0.489v3v7  - 0.843v5v6  +  0.0432x9x10  - 0.0556x9xn  - 0.000786x,2,  - 15.57 


51 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


52 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center 
Fort  Belvoir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  School 
Monterey,  California 

3.  Johannes  O.  Royset,  PhD 
Department  of  Operations  Research 
Monterey,  California 

4.  R.  Kevin  Wood,  PhD 
Department  of  Operations  Research 
Monterey,  California 

5.  Habib  Giirkan  Ba§ova 
Kara  Kuvvetleri  Komutanhgi 
Ankara,  TURKEY 

6.  Kara  Harp  Okulu  Kutuphanesi 
Bakanhklar  -  06100 
Ankara,  TURKEY 

7.  Savunma  Bilimleri  Enstitiisu 
Kara  Harp  Okulu 
Bakanhklar,  Ankara,  TURKEY 

8.  METU  Library 

Ortadogu  Teknik  Universitesi 
Inonu  Blv. 

Ankara,  TURKEY 

9.  Bilkent  University  Library 
Bilkent  Universitesi 
Ankara,  TURKEY 

10.  Ekrem  Erdem,  PhD 

Iktisadi  ve  Idari  Bilimler  Fakultesi 
Melikgazi,  Kayseri,  TURKEY 


53 


